import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import rasterio as rio
from rasterio.plot import show
from rasterio.plot import show_hist
import os
import fiona

!pwd

/c/Users/Bharath/anaconda3/envs/SG_Land_Reclamation_analysis

!ls

Approx_SG.kml
debug.log
Displaying images.ipynb
Extracting_SG_from_Landsat8_images.ipynb
Generated DataFrames
HDB data
Landsat_Image_processing_final_experiments.ipynb
Landsat8_image_processing_initial_experimentation.ipynb
Population_data_World_Bank.ipynb
Population_density_and_Housing_density_analysis.ipynb
Predicting the future of SG.ipynb
SG_cropped_Landsat8_images
SG_Landsat8_images
SG_shapefiles

os.chdir('SG_Landsat8_images/')

!pwd

/c/Users/Bharath/anaconda3/envs/SG_Land_Reclamation_analysis/SG_Landsat8_images

!ls

LC08_L1TP_125059_20130627_20170503_01_T1.tar
LC08_L1TP_125059_20140427_20170423_01_T1.tar
LC08_L1TP_125059_20141020_20170418_01_T1.tar
LC08_L1TP_125059_20150601_20170408_01_T1.tar
LC08_L1TP_125059_20160923_20170320_01_T1.tar
LC08_L1TP_125059_20170419_20170501_01_T1.tar
LC08_L1TP_125059_20180929_20181009_01_T1.tar
LC08_L1TP_125059_20190916_20190925_01_T1.tar
LC08_L1TP_125059_20200310_20200314_01_T1.tar

os.chdir('LC08_L1TP_125059_20200310_20200314_01_T1.tar/')
!ls

Bands_2020 = {}

for i in range(1,12):
    var = 'imgB'+str(i)
    
    img_file = 'LC08_L1TP_125059_20200310_20200314_01_T1_B' + str(i) + '.TIF'
    Bands_2020[var] = rio.open(img_file, 'r')

#print(Bands_2020)

os.chdir('..')

LC08_L1TP_125059_20200310_20200314_01_T1_ANG.txt
LC08_L1TP_125059_20200310_20200314_01_T1_B1.TIF
LC08_L1TP_125059_20200310_20200314_01_T1_B10.TIF
LC08_L1TP_125059_20200310_20200314_01_T1_B11.TIF
LC08_L1TP_125059_20200310_20200314_01_T1_B2.TIF
LC08_L1TP_125059_20200310_20200314_01_T1_B3.TIF
LC08_L1TP_125059_20200310_20200314_01_T1_B4.TIF
LC08_L1TP_125059_20200310_20200314_01_T1_B5.TIF
LC08_L1TP_125059_20200310_20200314_01_T1_B6.TIF
LC08_L1TP_125059_20200310_20200314_01_T1_B7.TIF
LC08_L1TP_125059_20200310_20200314_01_T1_B8.TIF
LC08_L1TP_125059_20200310_20200314_01_T1_B9.TIF
LC08_L1TP_125059_20200310_20200314_01_T1_BQA.TIF
LC08_L1TP_125059_20200310_20200314_01_T1_MTL.txt

os.chdir('LC08_L1TP_125059_20190916_20190925_01_T1.tar')
!ls

Bands_2019 = {}

for i in range(1,12):
    var = 'imgB'+str(i)
    
    img_file = 'LC08_L1TP_125059_20190916_20190925_01_T1_B' + str(i) + '.TIF'
    Bands_2019[var] = rio.open(img_file, 'r')

#print(Bands_2019)

os.chdir('..')

LC08_L1TP_125059_20190916_20190925_01_T1_ANG.txt
LC08_L1TP_125059_20190916_20190925_01_T1_B1.TIF
LC08_L1TP_125059_20190916_20190925_01_T1_B10.TIF
LC08_L1TP_125059_20190916_20190925_01_T1_B11.TIF
LC08_L1TP_125059_20190916_20190925_01_T1_B2.TIF
LC08_L1TP_125059_20190916_20190925_01_T1_B3.TIF
LC08_L1TP_125059_20190916_20190925_01_T1_B4.TIF
LC08_L1TP_125059_20190916_20190925_01_T1_B5.TIF
LC08_L1TP_125059_20190916_20190925_01_T1_B6.TIF
LC08_L1TP_125059_20190916_20190925_01_T1_B7.TIF
LC08_L1TP_125059_20190916_20190925_01_T1_B8.TIF
LC08_L1TP_125059_20190916_20190925_01_T1_B9.TIF
LC08_L1TP_125059_20190916_20190925_01_T1_BQA.TIF
LC08_L1TP_125059_20190916_20190925_01_T1_MTL.txt

os.chdir('LC08_L1TP_125059_20180929_20181009_01_T1.tar')
!ls

Bands_2018 = {}

for i in range(1,12):
    var = 'imgB'+str(i)
    
    img_file = 'LC08_L1TP_125059_20180929_20181009_01_T1_B' + str(i) + '.TIF'
    Bands_2018[var] = rio.open(img_file, 'r')

#print(Bands_2018)

os.chdir('..')

LC08_L1TP_125059_20180929_20181009_01_T1_ANG.txt
LC08_L1TP_125059_20180929_20181009_01_T1_B1.TIF
LC08_L1TP_125059_20180929_20181009_01_T1_B10.TIF
LC08_L1TP_125059_20180929_20181009_01_T1_B11.TIF
LC08_L1TP_125059_20180929_20181009_01_T1_B2.TIF
LC08_L1TP_125059_20180929_20181009_01_T1_B3.TIF
LC08_L1TP_125059_20180929_20181009_01_T1_B4.TIF
LC08_L1TP_125059_20180929_20181009_01_T1_B5.TIF
LC08_L1TP_125059_20180929_20181009_01_T1_B6.TIF
LC08_L1TP_125059_20180929_20181009_01_T1_B7.TIF
LC08_L1TP_125059_20180929_20181009_01_T1_B8.TIF
LC08_L1TP_125059_20180929_20181009_01_T1_B9.TIF
LC08_L1TP_125059_20180929_20181009_01_T1_BQA.TIF
LC08_L1TP_125059_20180929_20181009_01_T1_MTL.txt

os.chdir('LC08_L1TP_125059_20170419_20170501_01_T1.tar')
!ls

Bands_2017 = {}

for i in range(1,12):
    var = 'imgB'+str(i)
    
    img_file = 'LC08_L1TP_125059_20170419_20170501_01_T1_B' + str(i) + '.TIF'
    Bands_2017[var] = rio.open(img_file, 'r')

#print(Bands_2017)

os.chdir('..')

LC08_L1TP_125059_20170419_20170501_01_T1_ANG.txt
LC08_L1TP_125059_20170419_20170501_01_T1_B1.TIF
LC08_L1TP_125059_20170419_20170501_01_T1_B10.TIF
LC08_L1TP_125059_20170419_20170501_01_T1_B11.TIF
LC08_L1TP_125059_20170419_20170501_01_T1_B2.TIF
LC08_L1TP_125059_20170419_20170501_01_T1_B3.TIF
LC08_L1TP_125059_20170419_20170501_01_T1_B4.TIF
LC08_L1TP_125059_20170419_20170501_01_T1_B5.TIF
LC08_L1TP_125059_20170419_20170501_01_T1_B6.TIF
LC08_L1TP_125059_20170419_20170501_01_T1_B7.TIF
LC08_L1TP_125059_20170419_20170501_01_T1_B8.TIF
LC08_L1TP_125059_20170419_20170501_01_T1_B9.TIF
LC08_L1TP_125059_20170419_20170501_01_T1_BQA.TIF
LC08_L1TP_125059_20170419_20170501_01_T1_MTL.txt

os.chdir('LC08_L1TP_125059_20160923_20170320_01_T1.tar')
!ls

Bands_2016 = {}

for i in range(1,12):
    var = 'imgB'+str(i)
    
    img_file = 'LC08_L1TP_125059_20160923_20170320_01_T1_B' + str(i) + '.TIF'
    Bands_2016[var] = rio.open(img_file, 'r')

#print(Bands_2016)

os.chdir('..')

LC08_L1TP_125059_20160923_20170320_01_T1_ANG.txt
LC08_L1TP_125059_20160923_20170320_01_T1_B1.TIF
LC08_L1TP_125059_20160923_20170320_01_T1_B10.TIF
LC08_L1TP_125059_20160923_20170320_01_T1_B11.TIF
LC08_L1TP_125059_20160923_20170320_01_T1_B2.TIF
LC08_L1TP_125059_20160923_20170320_01_T1_B3.TIF
LC08_L1TP_125059_20160923_20170320_01_T1_B4.TIF
LC08_L1TP_125059_20160923_20170320_01_T1_B5.TIF
LC08_L1TP_125059_20160923_20170320_01_T1_B6.TIF
LC08_L1TP_125059_20160923_20170320_01_T1_B7.TIF
LC08_L1TP_125059_20160923_20170320_01_T1_B8.TIF
LC08_L1TP_125059_20160923_20170320_01_T1_B9.TIF
LC08_L1TP_125059_20160923_20170320_01_T1_BQA.TIF
LC08_L1TP_125059_20160923_20170320_01_T1_MTL.txt

!pwd

/c/Users/Bharath/anaconda3/envs/SG_Land_Reclamation_analysis/SG_Landsat8_images

Function to display Bands 1 - 9 for each Landsat Image

def show_bands(dictionary):
    
    Bands = dictionary
    fig, ax = plt.subplots(3,3, figsize = (15,15))

    num = 0
    for i in range(3):
        for j in range(3):
            num += 1
            key = 'imgB' + str(num)
        
            ax[i,j].imshow(Bands[key].read(1))
            ax[i,j].set_xlabel(key)
        

Function to display Histogram of Bands 1 - 9 for each Landsat image

def show_histogram(dictionary):
    Bands = dictionary
    fig, ax = plt.subplots(3,3, figsize = (15,15))
    
    num = 0
    for i in range(3):
        for j in range(3):
            num += 1
            key = 'imgB' + str(num)
        
            show_hist(Bands[key], bins = 50, ax = ax[i,j])
            #ax[i,j].set_xlabel(key)
    

!pwd

/c/Users/Bharath/anaconda3/envs/SG_Land_Reclamation_analysis/SG_Landsat8_images

os.chdir('..')

!pwd

/c/Users/Bharath/anaconda3/envs/SG_Land_Reclamation_analysis

show_bands(Bands_2020)

show_histogram(Bands_2020)

show_bands(Bands_2019)

show_histogram(Bands_2019)

show_bands(Bands_2018)

show_histogram(Bands_2018)

show_bands(Bands_2017)

show_bands(Bands_2016)

fig, ax = plt.subplots(1,2, figsize = (15,15))

ax[0].imshow(Bands_2016['imgB5'].read(1))
ax[1].imshow(Bands_2016['imgB1'].read(1))

<matplotlib.image.AxesImage at 0x2291075ed00>

plt.figure(figsize = (20,10))
plt.imshow((Bands_2016['imgB5'].read(1) - Bands_2016['imgB1'].read(1))[3700:4700, 1500:3700])

<matplotlib.image.AxesImage at 0x2290ef16400>

plt.figure(figsize = (20,10))
plt.imshow((Bands_2016['imgB5'].read(1) - Bands_2016['imgB4'].read(1))[3700:4700, 1500:3700])

<matplotlib.image.AxesImage at 0x2290f12d5e0>

plt.figure(figsize = (20,10))
plt.imshow((Bands_2017['imgB5'].read(1) - Bands_2017['imgB4'].read(1))[3700:4700, 1500:3700])

<matplotlib.image.AxesImage at 0x2293e86c640>

plt.figure(figsize = (20,10))
plt.imshow((Bands_2018['imgB5'].read(1) - Bands_2018['imgB4'].read(1))[3700:4700, 1500:3700])

<matplotlib.image.AxesImage at 0x2291725fa60>

plt.figure(figsize = (20,10))
plt.imshow((Bands_2019['imgB5'].read(1) - Bands_2019['imgB4'].read(1))[3700:4700, 1500:3700])

<matplotlib.image.AxesImage at 0x2293e7d64c0>

(Bands_2019['imgB5'].read(1)).shape

(7741, 7591)

(Bands_2019['imgB9'].read(1)).shape

(7741, 7591)

plt.figure(figsize = (15,15))
plt.imshow(Bands_2019['imgB5'].read(1))

<matplotlib.image.AxesImage at 0x2290ef448e0>

plt.figure(figsize = (15,15))
plt.imshow(Bands_2019['imgB9'].read(1))

<matplotlib.image.AxesImage at 0x22906fdb820>

Function to minimize the effects of Cloud Cover

This is a major challenge for analyzing Landsat images for locations like Singapore that are very cloudy almost throughout the year.

def remove_clouds(dictionary):
    Bands = dictionary
    
    img_without_cloud = Bands['imgB5'].read(1) - Bands['imgB9'].read(1)
    
    fig, ax = plt.subplots(1,3, figsize = (15,15))
    
    ax[0].imshow(Bands['imgB5'].read(1))
    ax[0].set_xlabel('imgB5')
    ax[1].imshow(Bands['imgB9'].read(1))
    ax[1].set_xlabel('imgB9')
    ax[2].imshow(img_without_cloud)
    ax[2].set_xlabel('img_without_cloud')
    
    return img_without_cloud
    

wc_2020 = remove_clouds(Bands_2020)

 plt.imshow(wc_2020)

<matplotlib.image.AxesImage at 0x2293e608430>

plt.figure(figsize=(20,10))
plt.imshow(wc_2020[3700:4700, 1500:3700])

<matplotlib.image.AxesImage at 0x2290ef61490>

wc_2019 = remove_clouds(Bands_2019)

 plt.imshow(wc_2019)

<matplotlib.image.AxesImage at 0x229170b0a00>

plt.figure(figsize=(20,10))
plt.imshow(wc_2019[3700:4700, 1500:3700])

<matplotlib.image.AxesImage at 0x22916ff20d0>

wc_2018 = remove_clouds(Bands_2018)

plt.figure(figsize=(20,10))
plt.imshow(wc_2018[3700:4700, 1500:3700])

<matplotlib.image.AxesImage at 0x22906fbce80>

wc_2017 = remove_clouds(Bands_2017)

plt.figure(figsize=(20,10))
plt.imshow(wc_2017[3700:4700, 1500:3700])

<matplotlib.image.AxesImage at 0x2290ef45a00>

wc_2016 = remove_clouds(Bands_2016)

plt.figure(figsize=(20,10))
plt.imshow(wc_2016[3700:4700, 1500:3700])

<matplotlib.image.AxesImage at 0x2290f30c190>

Function to find Normalized difference

1. NDWI : Normalized Difference Water Index (Uses Landsat8 Band 3 and Band 6 )
2. NDVI : Normalized Difference Vegetation Index (Uses Landsat8 Band 4 and Band 5 )

def NormalizedDiff(band1, band2, band9):
    
    b1_wo_clouds = band1.read(1)
    b2_wo_clouds = band2.read(1)
    
    #b1_wo_clouds = band1.read(1) - band9.read(1)
    #b2_wo_clouds = band2.read(1) - band9.read(1)
    
    b1 = b1_wo_clouds.astype('float')
    b2 = b2_wo_clouds.astype('float')
    
    ND = np.where(b1+b2==0., -1.0, (b1 - b2)/(b1 + b2))
        
    return ND
    

NDWI_2020 = NormalizedDiff(Bands_2020['imgB3'],Bands_2020['imgB6'], Bands_2020['imgB9'] )

<ipython-input-57-ba7b939c5a5e>:12: RuntimeWarning: invalid value encountered in true_divide
  ND = np.where(b1+b2==0., -1.0, (b1 - b2)/(b1 + b2))

plt.imshow(NDWI_2020)

<matplotlib.image.AxesImage at 0x229170d9ee0>

plt.figure(figsize=(20,10))
plt.imshow(NDWI_2020[3700:4800, 1500:3700])

<matplotlib.image.AxesImage at 0x22916e7b070>

NDWI_2019 = NormalizedDiff(Bands_2019['imgB3'],Bands_2019['imgB6'], Bands_2019['imgB9'])

<ipython-input-57-ba7b939c5a5e>:12: RuntimeWarning: invalid value encountered in true_divide
  ND = np.where(b1+b2==0., -1.0, (b1 - b2)/(b1 + b2))

plt.imshow(NDWI_2019)

<matplotlib.image.AxesImage at 0x2290f0955e0>

plt.figure(figsize=(20,10))
plt.imshow(NDWI_2019[3700:4700, 1500:3700])

<matplotlib.image.AxesImage at 0x229106ce160>